#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <complex.h>
#include <time.h>
#include <string.h>
#include "algo.h"

#define HALLOCSIZE 150000000000
void * hallocMem;
long int byteCounter = 0;
unsigned long long int mult=0, soma=0;
void hallocInit (void)
{
	hallocMem = (void *) malloc(HALLOCSIZE);
}

void hallocFree (void)
{
	free(hallocMem);
}

void * halloc (long int size)
{
	void * resp;

	resp = ((char *)hallocMem + byteCounter);
	byteCounter = (byteCounter + size) % (HALLOCSIZE - size);

	return resp;
}

void * challoc (long int size, int elementSize)
{
	void * resp;

	resp = ((char *)hallocMem + byteCounter);
	byteCounter = (byteCounter + size * elementSize) % (HALLOCSIZE - size * elementSize);

	memset(resp, 0, size * elementSize);

	return resp;
}

int nextPowerOfTwo (int n)
{
	int i = 0, original = n;
	int aux, lowerPower;

	while (n > 0)
	{
		n >>= 1;
		i++;
	};

	if (i == 0)
		return 1;

	aux = 2 << (i - 1);
	if (original % 2)
		return aux;
	else
	{
		lowerPower = nextPowerOfTwo(original - 1);
		if (lowerPower < aux)
			return lowerPower;
		else
			return aux;
	}
}

double * paddedVector (double *p, int originalSize)
{
	int i;
	int nSize;
	double *v;

	nSize = nextPowerOfTwo(originalSize);
	v = challoc(nSize, sizeof(double));

	for (i = 0; i < originalSize; i++)
		v[i] = p[i];

	return v;
}

Polynomial createPolynomial (int size)
{
	Polynomial new;

	new.size = size;

	new.p = halloc(size * sizeof(double));

	return new;
}

void freePolynomial (Polynomial p)
{
	free(p.p);
}


// retorna o máximo inteiro entre dois inteiros sem sinal
unsigned long int max( unsigned long int nA, unsigned long int nB){
	return nA > nB ? nA : nB;
}

// verifica se os polinomios A e B sao iguais testando-os em um numero qualquer
unsigned long int PolinomiosIguais(unsigned long int n, double A[], char *NA, double B[], char *NB)
{
	int r=1;
	for (unsigned long int i=0;i<n;i++) if (fabs(A[i]-B[i])>EPS) {
		printf("%10s: %lf.X^%lu\n",NA,A[i],i);
		printf("%10s: %lf.X^%lu\n\n",NB,B[i],i);
		r=0;
	}
	return(r);
}

/* Chamada do método direto e retorna o tempo gasto
 *
 */
unsigned long QuadCaller(unsigned long int nA, double A[],unsigned long int nB, double B[],double C[])
{
	//nA e nB são o números de termos dos respectivos polinômios associados
	unsigned long tinicio;
	int n;
	n = max(nA, nB);
	/* tratamento caso os polinomios não tenham a mesma ordem*/
	for (unsigned long int i=nA;i<n  ;i++) A[i] = 0.0; // completa com 0 posições restantes de A
	for (unsigned long int i=nB;i<n  ;i++) B[i] = 0.0; // completa com 0 posições restantes de B
	/////////////////////////////////////////////////////////////////////////////
	mult = 0;
	soma = 0;
	tinicio=clock();
	multpoli_quad(n, A, B, C);
	return(clock()-tinicio);
	/////////////////////////////////////////////////////////////////////////////
}



/* Chamada do método Karatsuba e retorna o tempo gasto
 *
 */
unsigned long KaraCaller(unsigned long int nA, double A[],unsigned long int nB, double B[], double C[])
{
	//nA e nB são o números de termos dos respectivos polinômios associados
	unsigned long tinicio;
	int n;
	n = max(nA, nB);
	/* tratamento caso os polinomios não tenham a mesma ordem*/
	for (unsigned long int i=nA;i<n  ;i++) A[i] = 0.0; // completa com 0 posições restantes de A
	for (unsigned long int i=nB;i<n  ;i++) B[i] = 0.0; // completa com 0 posições restantes de B
	/////////////////////////////////////////////////////////////////////////////
	mult = 0;
	soma = 0;
	tinicio=clock();
	multpoli_kara(n,A,B,C);
	return(clock()-tinicio);
	/////////////////////////////////////////////////////////////////////////////
}
/* Chamada do método FFT e retorna o tempo gasto
 *
 */
unsigned long FFTCaller(unsigned long int nA, double A[], unsigned long int nB, double B[], double C[])
{
	//nA e nB são o números de termos dos respectivos polinômios associados
	unsigned long tinicio;
	int n;
	n = max(nA, nB);
	/* tratamento caso os polinomios não tenham a mesma ordem*/
	for (unsigned long int i=nA;i<n  ;i++) A[i] = 0.0; // completa com 0 posições restantes de A
	for (unsigned long int i=nB;i<n  ;i++) B[i] = 0.0; // completa com 0 posições restantes de B
	/////////////////////////////////////////////////////////////////////////////
	mult = 0;
	soma = 0;
	tinicio=clock();
	multpoli_fft(n, A, B, C);
	return (clock()-tinicio);
}


int TestaPorArquivo(char *filename)
{
	unsigned long int np,nA,nB, i;
	unsigned long long int multquad=0, multkara=0, multfft=0, somaquad=0, somakara=0, somafft=0;
	double *A, *B, *C1, *C2, *C3;
	double comp1, comp2, comp3;
	FILE *ifp, *ofp;
	unsigned long long int  tQuad=0,tKara=0, tFFT=0, tm_Quad=0, tm_Kara=0, tm_FFT=0;
	ifp = fopen(filename,"r");
	if (ifp==NULL){
		printf("Erro para abrir arquivo %s\n.\n",filename);
		return(0);
	}
	fscanf(ifp,"%lu",&np);
	/*Laço que percorre e resolve uma multiplicação polinomial por vez e retorna os ciclos
	 * para cada um dos algoritmos expostos Ordinário, Karatsuba e FFT
	 */
	for (i=0;i<np;i++){
		// le o polinomio A com nA elementos e o polinomio B com nB elementos
		fscanf(ifp,"%lu",&nA);
		A = malloc(nA * sizeof(double));
		for (unsigned long int j=0;j<nA;j++) fscanf(ifp,"%lf",&A[j]);
		fscanf(ifp,"%lu",&nB);
		B = malloc(nB * sizeof(double));
		for (unsigned long int j=0;j<nB;j++) fscanf(ifp,"%lf",&B[j]);
		C1 = malloc((nA + nB) * sizeof(double));
		C2 = malloc((nA + nB) * sizeof(double));
		C3 = malloc((nA + nB) * sizeof(double));

		tQuad += QuadCaller(nA, A, nB, B, C1); 	// acumula os ciclos gastos pelo met. ordinario
		multquad += mult;
		somaquad += soma;
		tKara += KaraCaller(nA, A, nB,B, C2);  	// acumula os ciclos gastos pelo met. Karatsuba
		multkara += mult;
		somakara += soma;
		tFFT += FFTCaller(nA, A, nB, B, C3); 	// acumula os ciclos gastos pelo met. FFT
		multfft += mult;
		somafft += soma;

		tm_Quad=+tm_Quad+tQuad;
		tm_Kara=+tm_Quad+tKara;
		tm_FFT=+tm_FFT+tFFT;
	}
	tm_Quad=+tm_Quad/np;
	tm_Kara=+tm_Quad/np;
	tm_FFT=+tm_FFT/np;

	fclose(ifp);

	/*Retorna o tempo medio de execucao de um arquivo de testes em um arquivo de saida*/
	ofp = fopen("ra035501_136640.log", "a+");//cria arquivo se não existir e realiza um update nele a partir do fim do arquivo
	if(ofp == NULL){
		fprintf(stderr, "Nao foi possivel abrir o arquivo de saida!\n", "ra035501_136640.log");
		exit(1);
	}
	fprintf(ofp,"quad %s %llu %llu %llu\n", filename, multquad, somaquad, tQuad);
	fprintf(ofp,"kara %s %llu %llu %llu\n", filename, multkara, somakara, tKara);
	fprintf(ofp,"fft %s %llu %llu %llu\n", filename, multfft, somafft, tFFT);
	fclose(ofp);
	if(tQuad <= tKara && tQuad <= tFFT) return 0;
	if(tKara < tQuad && tKara < tFFT) return 1;
	if(tFFT < tQuad && tFFT>tQuad) return 2;
	else return -1;

}

int TestaAleatorio(unsigned long int n)
{
	double *A, *B, *C1, *C2, *C3;
	unsigned long long int multquad=0, multkara=0, multfft=0, somaquad=0, somakara=0, somafft=0;

	A = malloc(n * sizeof(double));
	B = malloc(n * sizeof(double));
	C1 = malloc(2 * n * sizeof(double));
	C2 = malloc(2 * n * sizeof(double));
	C3 = malloc(2 * n * sizeof(double));

	double comp1, comp2, comp3;
	unsigned long long int tQuad=0,tKara=0, tFFT=0;
	srand(clock()); // gera dois polinomios aleatorios de tam. n
	for (int i=0;i<100;i++) {// 100 pares de polinomios
		for (unsigned long int j=0;j<n;j++){// cada polinomio tem tamanho n
			A[j]=(double)((int)10*rand());
			B[j]=(double)((int)10*rand());
		}
		///////////////////////////////////////////////////////////////////
		tQuad += QuadCaller(n, A, n, B, C1);// resolve por met. DummyMult
		multquad += mult;
		somaquad += soma;
		tKara += KaraCaller(n, A, n, B, C2);  // resolve por met. Karatsuba
		multkara += mult;
		somakara += soma;
		tFFT += FFTCaller(n, A, n, B, C3); //resolve por met. FFT
		multfft += mult;
		somafft += soma;
		///////////////////////////////////////////////////////////////////
		if (!PolinomiosIguais(n+n-1,C1,"Ordinario",C2,"Karatsuba")) {
			printf("\n\nIter. %d: Polinomios sao diferentes.\n",i+1);
			printf("Polinomios C1 e C2 sao diferentes.\n"); break;}
		if (!PolinomiosIguais(n+n-1,C1,"Ordinario",C3,"FFT")) {
			printf("\n\nIter. %d: Polinomios sao diferentes.\n",i+1);
			printf("Polinomios C1 e C3 sao diferentes.\n"); break;}
		if (!PolinomiosIguais(n+n-1,C2,"Karatsuba",C3,"FFT")) {
			printf("\n\nIter. %d: Polinomios sao diferentes.\n",i+1);
			printf("Polinomios C2 e C3 sao diferentes.\n"); break;}
		printf("\nIt.%4u: Tempo Ordinario = %llu, Tempo Karatsuba = %llu, Tempo FFT = %llu \nOrdinario/Karatsuba = %5.3lf, Ordinario/FFT = %5.3lf, Karatsuba/FFT = %5.3lf\n", i+1, tQuad, tKara, tFFT, (double)tQuad/(double)tKara,(double)tQuad/(double)tFFT, (double)tKara/(double)tFFT);
		printf("Quad Mult = %llu Somas = %llu\n", multquad, somaquad);
		printf("Kara Mult = %llu Somas = %llu\n", multkara, somakara);
		printf("FFT Mult = %llu Somas = %llu\n", multfft, somafft);
	}
	/*Retorna o algoritmo mais rápido*/
	///Compara Ordinario contra Karatsuba
	comp1 = (double)tQuad/(double)tKara;
	///Compara Karatsuba contra FFT
	comp2 = (double)tKara/(double)tFFT;
	///Compara Ordinario contra FFT
	comp3 = (double)tQuad/(double)tFFT;

	if(tQuad <= tKara && tQuad <= tFFT) return 0;
	if(tKara < tQuad && tKara < tFFT) return 1;
	if(tFFT < tQuad && tFFT< tQuad) return 2;
	else return -1;//convenciona erro
}


int main(int argc,char **argv)
{
	hallocInit();
	unsigned long int n; int resultado;
	if ((argc>2)||(argc<1)) {printf("Numero de parametros errado");return(1);}
	if (argc==1) {
		printf("Ha' duas alternativas de uso:\n\n"
				"                   %s <arquivo>\n\n"
				"                   %s <tamanho_polinomios>\n\n"
				"Na primeira, <arquivo> e' nome do arquivo onde os dados serao lidos e testados.\n"
				"Na segunda, <tamanho_polinomios> e' o numero de coeficientes dos\n"
				"polinomios a multiplicar (serao executados 100 testes aleatorios).\n",
				argv[0],argv[0]);
		return(1);
	}
	else {
		n = atoi(argv[1]);
		if (n==0){
			resultado = TestaPorArquivo(argv[1]);
			hallocFree();
		}
		else{
			resultado = TestaAleatorio(n);
			hallocFree();
		}
	}

	if (resultado == 0) printf("\nMultiplicacao direta foi mais rapida, na media.\n\n");
	if (resultado == 1) printf("\nMultiplicacao por Karatsuba foi mais rapida, na media.\n\n");
	if (resultado == 2) printf("\nMultiplicacao por FFT foi mais rapida, na media.\n\n");
	if (resultado == -1) printf("\nERRO!\n\n");

//	hallocFree();

	return(0);

}

complex * doubleToComplex (double *p, int size, int biggerSize)
{
	complex *newV = challoc(biggerSize, sizeof(complex));
	int i;

	for (i = 0; i < size; i++)
		newV[i] = p[i];

	return newV;
}

double * complexToDouble (complex *p, int size)
{
	double *newV = halloc(size * sizeof(double));
	int i;

	for (i = 0; i < size; i++)
		newV[i] = creal(p[i]);

	return newV;
}

/**
 * Dado um vetor de complexos v de tamanho N, retorna uma cópia sua.
 */
complex* copyCVec (complex * v, int N)
{
	complex *newV = halloc(N * sizeof(complex));
	int i;

	for (i = 0; i < N; i++)
		newV[i] = v[i];

	return newV;
}

/**
 * Dado um vetor v de tamanho N, retorna uma cópia sua.
 */
double * copyVec (double * v, int N)
{
	double *newV = halloc(N * sizeof(double));

	memcpy(newV, v, N * sizeof(double));

	return newV;
}

/**
 * Dada uma representação polinomial p de N termos, retorna a primeira
 * parcela da decomposição de p (envolvendo somente os termos de ordem
 * par).
 */
complex * computeP1 (complex * p, int N)
{
	complex * p1 = challoc(N / 2, sizeof(complex));
	int i;

	for (i = 0; i < N / 2; i++)
		p1[i] = p[2 * i];

	return p1;
}

/**
 * Dada uma representação polinomial p de N termos, retorna a segunda
 * parcela da decomposição de p (envolvendo somente os termos de ordem
 * ímpar).
 */
complex * computeP2 (complex * p, int N)
{
	complex * p2 = challoc(N / 2, sizeof(complex));
	int i;

	for (i = 0; i < N / 2; i++)
		p2[i] = p[2 * i + 1];

	return p2;
}

/**
 * Função que divide todos os N elementos de um vetor v por scale.
 * Deve ser usada pela função de FFT inversa.
 */
void vecScale (complex * v, int N, int scale)
{
	int i;

	for (i = 0; i < N; i++)
	{
		v[i] = v[i] / scale;
	}
}



/**
 * Função de operação de FFT sobre uma representação polinomial p
 * de N termos.
 */
complex * recDFT (complex * p, int N)
{
	complex WN = cexpf((2.0f / N) * M_PI * I);
	complex W = 1.0f;
	complex *p1, *p2;
	complex *d1, *d2, *d;
	int k;

	if (N == 1)
		return copyCVec(p, 1);

	p1 = computeP1(p, N);
	p2 = computeP2(p, N);
	d1 = recDFT(p1, N / 2);
	d2 = recDFT(p2, N / 2);

	d = challoc(N, sizeof(complex));

	for (k = 0; k < N / 2; k++)
	{
		d[k] = d1[k] + W * d2[k];
		d[k + N / 2] = d1[k] - W * d2[k];
		W = W * WN;
		soma = soma + 2;
		mult = mult + 2;
	}

	return d;
}

/**
 * Função de operação de FFT inversa sobre uma representação polinomial
 * p de N termos. O resultado é reescalado corretamente.
 */
complex * recIDFT (complex * p, int N)
{
	complex WN = cexp(-(2.0f / N) * M_PI * I);
	complex W = 1.0f;
	complex *p1, *p2;
	complex *d1, *d2, *d;
	int k;

	if (N == 1)
		return copyCVec(p, N);

	p1 = computeP1(p, N);
	p2 = computeP2(p, N);
	d1 = recIDFT(p1, N / 2);
	d2 = recIDFT(p2, N / 2);

	d = challoc(N, sizeof(complex));

	for (k = 0; k < N / 2; k++)
	{
		d[k] = d1[k] + W * d2[k];
		d[k + N / 2] = d1[k] - W * d2[k];
		W = W * WN;
		mult=mult+2;
		soma=soma+2;
	}
	return d;
}

/**
 * Wrapper de chamada da função FFT.
 */
complex * evaluateDFT (double * p, int N)
{
	complex *ccopy, *result;

	ccopy = doubleToComplex(p, N, 2 * nextPowerOfTwo(N));
	result = recDFT(ccopy, 2 * nextPowerOfTwo(N));

	return result;
}

/**
 * Função de multiplicação ponto a ponto.
 * Dadas duas representações polinomiais p1 e p2 de N termos
 * cada uma, retorna um novo polinômio cujos coeficientes correspondem
 * ao produto escalar entre os vetores de coeficientes de p1 e p2.
 */
complex * pointwiseMult (complex * p1, complex * p2, int N)
{
	complex * p = challoc(N, sizeof(complex));
	int i;

	for (i = 0; i < N; i++)
		p[i] = p1[i] * p2[i];
	/* Esta operação conta como multiplicação na computação total do FFT*/
	mult++;

	return p;
}

/**
 * Função de multplicação de polinômios segundo algoritmo clássico.
 * Dadas duas representações polinomiais p1 e p2 de n termos cada,
 * retorna o produto de tais polinômios.
 */
double * quad_algo (int n, double * p1, double * p2)
{
	double * p = challoc(2 * n - 1, sizeof(double));
	int i, j;

	for (i = 0; i < n; i++)
		for (j = 0; j < n; j++){
			p[i + j] += p1[i] * p2[j];
			mult++;
			soma = mult-1;
		}
	return p;
}

/**
 * Função de multplicação de polinômios segundo algoritmo FFT.
 * Dadas duas representações polinomiais p1 e p2 de n termos cada,
 * retorna o produto de tais polinômios.
 */
double * fft_algo (int n, double * p1, double * p2)
{
	complex *d1, *d2, *d, *p;
	int nSize;

	d1 = evaluateDFT(p1, n);
	d2 = evaluateDFT(p2, n);

	nSize = 2 * nextPowerOfTwo(n);

	d = pointwiseMult(d1, d2, nSize);

	p = recIDFT(d, nSize);

	vecScale(p, nSize, nSize);

	return complexToDouble(p, nSize);
}

/**
 * Função de determinação do ponto de split para o algoritmo Karatsuba.
 * Retorna um ponto mé
 */
int split_size (int N)
{
	return N / 2;
}

/**
 * Função de split de polinômios segundo método Karatsuba.
 */
struct SplitPolynomial split (int N, int split_point, struct Polynomial p)
{
	struct SplitPolynomial sp;
	double * aux;
	int i;

	sp.split = split_point;

	sp.high.p = halloc(N / 2 * sizeof(double));
	sp.high.size = split_point;
	sp.low.p  = halloc(N / 2 * sizeof(double));
	sp.low.size  = split_point;

	aux = p.p + N / 2;

	for (i = 0; i < N / 2; i++)
	{
		sp.low.p[i] = p.p[i];
		sp.high.p[i] = aux[i];
	}

	return sp;
}

/**
 * Função de subtração de polinômios.
 * Dados dois polinômios x e y, faz com que x contenha x - y.
 */
void sub_coeff (struct Polynomial x, struct Polynomial y)
{
	int i;

	if (x.size != y.size)
		printf("sub_coeff cannot handle different polynomial sizes.\n");

	for (i = 0; i < x.size; i++)
	{
		x.p[i] -= y.p[i];
	}
}

/**
 * Função de adição de polinômios.
 * Dados dois polinômios x e y, faz com que x contenha x + y.
 * Assume-se que os dois polinômios possuam mesmo número de termos.
 */
void add_coeff (struct Polynomial x, struct Polynomial y)
{
	int i;

	if (x.size != y.size)
		printf("add_coeff cannot handle different polynomial sizes.\n");

	for (i = 0; i < x.size; i++)
	{
		x.p[i] += y.p[i];
	}
}

Polynomial shiftPoly (Polynomial p, int amount)
{
	Polynomial new;
	int i;

	new.p = halloc((amount + p.size) * sizeof(double));
	new.size = p.size + amount;

	for (i = 0; i < amount; i++)
		new.p[i] = 0;

	for (i = amount; i < new.size; i++)
		new.p[i] = p.p[i];

	return new;
}

/**
 * Função de composição de polinômios a partir de sua representação
 * tripartite segundo relação p = z2*B^(2m) + z1*B^(m) + z0.
 */
struct Polynomial compose_poly (struct Z zees)
{
	struct Polynomial resp;
	int i, shift1, shift2;

	shift1 = zees.split;
	shift2 = 2 * zees.split;

	resp.size = zees.z2.size + 2 * zees.split;
	resp.p = challoc(resp.size, sizeof(double));

	for (i = 0; i < resp.size; i++)
	{
		if (i < zees.z0.size)
			resp.p[i] += zees.z0.p[i];

		if (shift1 <= i && i < shift1 + zees.z1.size)
			resp.p[i] += zees.z1.p[i - shift1];

		if (shift2 <= i && i < shift2 + zees.z2.size)
			resp.p[i] += zees.z2.p[i - shift2];
	}

	return resp;
}

/**
 * Função que calcula as três componentes z2, z1 e z0 da decomposição
 * polinomial empregada no algoritmo Karatsuba.
 */
struct Z computeZ (struct SplitPolynomial X, struct SplitPolynomial Y)
{
	struct Z resp;
	resp.z2.p = NULL;
	Polynomial aux;

	resp.z0 = karat_core(X.split, X.low, Y.low);
	resp.z2 = karat_core(X.split, X.high, Y.high);

	//X.high e Y.high passam a conter as respectivas somas
	add_coeff(X.high, X.low);
	add_coeff(Y.high, Y.low);

	resp.z1 = karat_core(X.split, X.high, Y.high);

	//Cria polinômio auxiliar que evita a sobrescrição de resp.z2 por add_coeff
	aux.size = resp.z2.size;
	aux.p = copyVec(resp.z2.p, aux.size);
	add_coeff(aux, resp.z0);

	sub_coeff(resp.z1, aux);

	resp.split = X.split;

	return resp;
}

/**
 * Função de liberação de memória de estrutura de struct Z.
 */
void freeZ (struct Z zees)
{
	freePolynomial(zees.z0);
	freePolynomial(zees.z1);
	freePolynomial(zees.z2);
}

/**
 * Função de multplicação de polinômios segundo algoritmo Karatsuba.
 * Dadas duas representações polinomiais p1 e p2 de n termos cada,
 * retorna o produto de tais polinômios.
 * Para n inferior a 10, o algoritmo clássico (quadrático) é empregado.
 */
struct Polynomial karat_core (int n, struct Polynomial x, struct Polynomial y)
{
	int split_point;
	struct SplitPolynomial X, Y;
	struct Polynomial resp;
	struct Z zees;

	if (n < 70)
	{
		resp.p = quad_algo(n, x.p, y.p);
		resp.size = 2 * n - 1;
		return resp;
	}

	split_point = split_size(n);

	X = split(n, split_point, x);		//Aloca memória:			N
	Y = split(n, split_point, y);		//Aloca memória:			N

	zees = computeZ(X, Y);				//Aloca memória: 			S

	resp = compose_poly(zees);			//Aloca memória:			S

	return resp;
}

double * karat_algo (int n, double * p1, double * p2)
{
	Polynomial P1, P2;
	int newSize;

	newSize = nextPowerOfTwo(n);

	P1.size = newSize;
	P2.size = newSize;

	if (newSize != n)
	{
		P1.p = paddedVector(p1, n);
		P2.p = paddedVector(p2, n);
	}
	else
	{
		P1.p = p1;
		P2.p = p2;
	}

	return karat_core(nextPowerOfTwo(n), P1, P2).p;
}

void multpoli_quad(int n, double A[], double B[], double C[])
{
	C=quad_algo(n, A, B);
}

void multpoli_kara(int n, double A[], double B[], double C[])
{
	C=karat_algo(n, A, B);
}

void multpoli_fft (int n, double A[], double B[], double C[])
{
	C=fft_algo(n, A, B);
}
